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ABSTRACT 

We investigate the formation of the first stars at the end of the cosmic dark ages 

■ with a suite of three-dimensional, moving mesh simulations that directly resolve the 

collapse of the gas beyond the formation of the first protostar at the centre of a dark 

Q . matter minihalo. The simulations cover more than 25 orders of magnitude in density 

^ ■ and have a maximum spatial resolution of 100 km, which extends well below the radius 

^ \ of individual protostars and captures their interaction with the surrounding gas. In 

C3 ■ analogy to previous studies that employed sink particles, we find that the Kepler ian 

disc around the primary protostar fragments into a number of secondary protostars, 
which is enabled by H2 collisional dissociation cooling and collision-induced emission. 
^ ■ The further evolution of the protostellar system is characterised by strong gravitational 

\ torques that transfer angular momentum between the secondary protostars formed in 

i/^ . the disc and the surrounding gas. This leads to the migration of about half of the 

10 ' secondary protostars to the centre of the cloud in a free-fall time, where they merge 

10 \ with the primary protostar and facilitate its growth to about five times the mass of 

^sq ■ the second most massive protostar. By the same token, a fraction of the protostars 

■ obtain angular momentum from other protostars via N-body interactions and migrate 
! to higher orbits. On average, only every third protostar survives until the end of the 

simulation. However, the number of protostars present at any given time increases 
monotonically, suggesting that the system will continue to grow beyond the limited 
period of time simulated here. 
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1 INTRODUCTION 

The formation of the first stars, also termed Popula- 
tion III (or Pop III), heralded an end to the physical 
simplicity of the Universe at the end of the cosmic dark 
ages (e.^. Barkana k, Loebl I2OOII: iBromm Larson| 2004| : 



Ciardi Ferraral l2005l : iGloverl l2005l : iBromm et all boool : 
Loebl I2OIQI ). Their radiation heated and ionised the inter- 



galac t ic medium (IGM; e . g. iBromm e t al. 2001 



'Schaerer 



I2OO2I : lAlvarez et aD I2OO6I : I Abel et al.lT2007. : .Johnson et al.. 
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20071 : lYoshida et all l2007l : IWhalen et all I2OO8I : iGreif et al.1 

200y), and the violent supernova explosions that marked 
their deaths enriched the primordial gas with the first heavy 
elements (e.g. iHegeret al.1 l2003l: lUmeda &: Nomotd l2003l : 



IWise Abe]li2008l - lGreif et al.ll20ld ). One of the main goals 
is therefore to derive their characteristics, such as their ini- 
tial mass function (IMF) and rotation rate, which govern 
their influence on the IGM and subsequent stellar genera- 
tions. 

With the exception of magnetic fields, whose strength 
are poorly understood, the initial conditions for this prob- 
lem are well-posed, since the fluctuation power spectrum 
and elemental composition of the gas after recombination 
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have been constrained by obs ervations of the cosm ic mi- 
crowave background (CMB; e.g. Komat su et aP bOOQV How- 
ever, the complex interplay of the relevant physical processes 
over many orders of magnitude in spatial and temporal scale 
have so far precluded a comprehensive answer. Nevertheless, 
pioneering one-dimensional calculations as well as three- 
dimensional numerical simulations have suggested that the 
first stars formed in dark matter (DM) minihaloes of mass 
Mvir ^ 10^-10^ M© at redshifts z > 20 due t o the formation 
of molecular hydrogen fHaiman et al." '1996; Teemarketal 
[l997: Omukai & Nishi 1998; Abel et al. 2002; Bromm et al 
[2002). The virialised gas cools to a minimum temperature of 
about 200 K, which is set by the internal rovibrational tran- 
sitions of H2. The gas cloud then becomes Jeans-unstable 
and undergoes runaway collapse, resulting i n the formation 
of a p rotostar at the centre of the halo ( Yosh ida et al.ll2006l , 
I2OO8I ). In the absence of fragmentation, an extrapolation of 
the high accretion rate within the Kelvin-Helmholtz time 
yields a characteristic mass of M* > 100 Mq once the star 
reaches the main sequence, two orders of magnitude higher 
than for present-day stars (|Kroupal l200l [Chabrier 2003). 

While providing a plausible estimate of the expected 
mass scale, this argument glosses over details of the collapse 
and accretion process. In particular, the Courant constraint 
limited the simulations to a narrow timespan centred on the 
formation of the first protostellar core. Fragmentation there- 
fore only appeared if the collaps e of a secondary clump was 
synchronised with the first fe.g. iTurk et al.|[2009l l. Further- 
more, the formation and possible fragmentation of a circum- 
stellar disc could not be followed. Recent studies have cir- 
cumvented this limitation by employing so- called sink parti- 
cles that represent growing protostars (e.g. iBate et aL 1 19951: 
Krumholz et al.ll2004l : Ijappsen et al]l2005l : iFederrath et all 



201QI ). Irrespective of the scale on which the sink particles 



were inserted, these studies found that the protostellar disc 
fragmented strongly into a group of protostars, due to the 
high accretion rate from the envelope onto the disc, and the 
very efficient cooling of the disc by H2 lines and collision- 
induced emission (Clark et al. 2008, 2011a,b; Stacy et al. 
2010, Greif et al. 2011). This led to the inference that Pop III 
stars may have formed with a range of masses, and that the 
final mass function depends on the efficiency of fragmenta- 
tion and merging. 

Unfortunately, the utilization of sink particles also 
comes at a cost. In particular, the boundary conditions im- 
posed on the gas near the accretion radius are not a self- 
consistent extension of the underlying hydrodynamic equa- 
tions. With the exception of supersonic inflows, where the 
trajectories become ballistic, the boundary between the sink 
particles and the surrounding gas may be subject to numer- 
ical artefacts. The complicated torques near the accretion 
radius are not captured accurately, which might artiflcially 
enhance or prevent fragmentation. The loss of resolution on 
the scale of the accretion radius also sets a minimum scale 
on which fragmentation may occur. Finally, close encoun- 
ters between sink particles may result in dynamical ejections 
(e.g. Greif 2011; Smith et al. 2011a,b), although the orbital 
energy may in fact be dissipated through unresolved torques. 
An accurate account of fragmentation and protostellar inter- 
actions on scales comparable to the physical sizes of the pro- 
tostars therefore requires a more self- consistent treatment. 

In the present study, we address these issues by per- 



forming three-dimensional, moving-mesh simulations of the 
formation and evolution of primordial prot ostellar systems 
in minihaloes. Similar to I Greif et al.l (|201lh , we start from 
cosmological initial conditions and follow the collapse of the 
gas up to the formation of the first protostar. Once the col- 
lapse stalls, we deviate from previous work and do not insert 
sink particles. Instead, we self-consistently evolve the chem- 
istry and hydrodynamics of the gas beyond the 'Courant 
myopia' of the first protostar to capture the formation and 
fragmentation of the circumstellar disc, as well as the evolu- 
tion of the nascent protostellar system for many dynamical 
times. The maximum spatial resolution of 100 km allows us 
to resolve the immediate vicinity of the protostars and their 
complicated interaction with the surrounding gas cloud. 

The structure of our work is as follows. In Section 2, we 
describe the numerical setup of the simulations, the chem- 
istry and cooling model, and our method for finding pro- 
tostars and determining their properties. In Section 3, we 
present the simulations and analyse the initial collapse of 
the minihaloes, the formation and fragmentation of the disc, 
the evolution of the protostellar system as a whole, and the 
merger and survival rates of the secondary protostars formed 
in the disc. Finally, in Section 4 we summarise our results 
and draw conclusions. All distances are quoted in proper 
units, unless noted otherwise. 



2 NUMERICAL METHODOLOGY 

The initial conditions of the simulations are taken from 
iGreif et al.1 (|201lh . which are cosmological simulations that 
employ multiple levels of refinement to follow the collapse 
of the central gas clouds in minihaloes. We briefiy review 
these simulations and describe the additions made to the 
code to follow the collapse of the gas well into the optically 
thick regime at nn ^ 10^^ cm~^. All of our simulations were 
performed with th e moving-mesh code AREPO, which is de- 
scribed in detail in ISpringell (|201Q| ). 



2.1 Dark matter simulations 

We first use a series of DM simulations to determine the 
locations of the minihaloes that are later fiagged for re- 
finement. The employed boxes have a side length of 250 
and 500 kpc (comoving), and use 128^ and 256^ particles 
of mass Mdm — 272 M© . We use a comoving gravitational 
softening length of 68 pc and imprint an initial fiuctuation 
power spectrum for a A cold dark matter (ACDM) cosmol- 
ogy at a starting redshift of 2; = 99 with matter density 
= 1 - = 0.27, baryon density Qh = 0.046, Hubble 
parameter h = ifo/100 kms~^ Mpc~^ = 0.71 (where Hq is 
the present Hu bble expansion rate) , and a spectral index 
Us = 0.96 fe.g. iKomatsu et "aP l2009l V The linearly extrap- 
olated present-day normalisation ag is varied with respect 
to the measured value of 0.81 to represent rare overdense 
regions of the universe that cannot be captured with the 
limited box sizes employed here. The DM simulations are 
evolved until the first minihalo with a virial mass exceeding 
5 X 10^ M0 collapses. The parameters of the simulations as 
well as the properties of the minihaloes are summarised in 
Table 1. 
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Table 1. Simulation Parameters 



Simulation 


Size 


Particles 




^dm,ref 


Mgas 






Tvir 


A 


^^coll 




(kpc) 




(Mo) 


(Mo) 


(Mo) 




(Mo) 


(pc) 






MHl 


500 


256^ 


272 


3.53 


0.72 


0.9 


3.0 X 10^ 


110 


0.055 


19.5 


MH2 


250 


1283 


272 


3.53 


0.72 


1.2 


2.3 X 10^ 


94 


0.073 


20.9 


MH3 


500 


256^ 


272 


3.53 


0.72 


1.1 


3.1 X 10^ 


97 


0.044 


22.6 


MH4 


500 


256^ 


272 


3.53 


0.72 


1.3 


1.8 X 10^ 


58 


0.038 


31.7 



The comoving box sizes, initial DM particle numbers, DM masses, refined DM masses, gas masses, and nor- 
malisations used in the cosmological simulations, followed by the virial masses, virial radii, spin parameters, 
and collapse redshift s of the first minihaloes that form. The ha l o pro perties agr e e well with those found in 
previous studies fe.g. iMachacek et al.ir200ll : lYoshida et al.llioosi : iGao et al.ll2007l : lO'Shea fc Normanll2007l V 



2.2 Refined simulations 



2.4 (De-) refinement strategy 



After the location of a minihalo has been determined, we 
construct a high-resolution study of the halo by first identi- 
fying all particles within approximately two times the virial 
radius. These particles are then traced back to their initial 
conditions, yielding the Lagrangian region that formed the 
halo. To construct new initial conditions, each particle in 
this region is replaced by 64 DM particles and 64 mesh- 
generating points that generate the space-filling cells that 
are used for the hydrodynamic calculations. 

Cells and DM particles outside of the high-resolution 
region are replaced by ever higher-mass particles and cells 
with increasing distance from the target halo, such that the 
resolution coarsens significantly far away from the target 
halo, yet the gravitational tidal field that influences its for- 
mation is still followed with high accuracy. This step reduces 
the total number of particles and cells in the refined initial 
conditions to :^ 2 x 10^. The simulation box is then centred 
on the target halo and reinitialised at z = 99 based on the 
original realisation of the density field, but augmented with 
additional small-scale power in the high-resolution region 
that can now be represented. The initial DM particle and 
cell masses in the high-resolution region before any further 
run-time refinement are Mdm,ref = (1 — nb/^m)Mdm/64 ~ 
3.53 Mo and Mgas (nb/^^m)Mdm/64 :^ 0.72 Mo, respec- 
tively. We use a comoving gravitational softening length of 
17 pc for the refined DM component. 



2.3 Extraction procedure 

The refined simulations are run until the first cell exceeds a 
density of nn = 10^ cm"*^, after which we extract the central 
1 pc and reinitialise the simulations with reflective bound- 
ary conditions, which are better suited to deal with highly 
asymmetric boundaries. We further remove all DM particles, 
which has almost no effect on the evolution of the central gas 
cloud, since the gas is already well decoupled from the DM 
component (|Greif et al.1 l201lh . These initial resimulations 
are then evolved until the density exceeds nn = 10^^ cm""^, 
after which we again extract the central 2000 AU and employ 
reflective boundary conditions. The sound crossing time of 
the refined box exceeds the simulation time by at least two 
orders of magnitude, such that perturbations from the edges 
of the box do not affect our results. 



The collapse of the gas to ever smaller scales leads to a 
continuous decrease of the Jeans mass, which must be ac- 
companied by a continuous refinement of the gas to ensure 
that the Truelove criterion is fulfilled (Truelove et al. 1991). 
W e use a mod i fied ve rsion of the Jeans refinement described 
m iGreif et al. 1 (|201lh . which refines a cell once its density 
exceeds nn = lcm~^ within the central 200 pc of the halo, 
and if in addition the Jeans length falls below the target 
number of cells per Jeans length A/j, which is given by 



log (nH/nH,a 
log (nH,b/nH,< 



10^cm-^ nH,b = 10^^cm-^ Aj,, 



(1) 



128, 



Aj,a, and for 



where nn,; 

and Aj,b = 32. For nu < nH,a, Aj 
nn > nH,b, Aj = Aj,b- The utilisation of 128 cells per Jeans 
length at low densities ensures that the turbulence gener- 
ated during_th^_Jnitm adequately 
(e.g. iFederrath et"al]l201ll : iTurk et al.ll2oTi ) , while the lower 
resolution of 32 cells per Jeans length at high densites en- 
sures that the simulations remain computationally feasible. 
We have found that even a slight increase of Aj,b to 64 results 
in an increase of the necessary computing time by nearly an 
order of magnitude, which cannot be easily compensated by 
using a larger number of processors. Each simulation was 
run for about three months on 32 processors, and employed 
several million internal timesteps to evolve the simulations 
for about 10 yr. 

The cells that are flagged for reflnement are split into 
two cells by replacing the mesh generating point with a close 
pair of points. The conserved fluid quantities in the original 
cell are distributed onto the two halves in proportion to 
their volume fractions. In the subsequent evolution, the pair 
of mesh-generating points separates, leading to a smooth 
adjustment of the local mesh. We also use the dereflnement 
procedure of AREPO to save computational time in regions 
where the collapse has stalled. If the size of a cell exceeds the 
Jeans length by a factor of two, and if in addition it is part 
of a smooth region of the flow (measured by requiring that 
the density, velocity, pressure, and temperature gradients 
across a cell are small), we remove the cell's mesh-generating 
point and distribute its content onto the neighboring cells 
in a conservative fashion. The relative fraction given to each 
of the neighboring cells is determined by the fraction of the 
volume that they claim from the removed cell, as determined 
by the Voronoi tessellation used in AREPO. 
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Figure 1. From top left to bottom right: density of hydrogen nuclei, enclosed gas mass, temperature, H2 fraction, e~ fraction, radial 
velocity, sound speed, radial Mach number, rotation velocity, ratio of rotation to Keplerian velocity, angle between the vector of the 
rotation velocity at a given radial bin and the innermost radial bin, and turbulent Mach number versus distance to the densest cell in the 
simulation. The various line styles denote the minihaloes according to the legend. The above quantities are calculated when the density 
first exceeds nn = 10^^ cm~^ and a shock forms at the surface of the primary protostar, which is located at about 0.01 AU. The plot 
covers all relevant length scales and shows the collapse of the gas from cosmological to protostellar scales (see Section 3.1 for a detailed 
analysis) . 
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2.5 Chemistry and cooling 

We have substantiah y revis ed the chemical network em- 
ployed in Greif et alJ (|201lh to include a non-equilibrium 
solver at low densities and an equilibrium solver at high 
densities. The non-equilibrium solver uses the publicly avail- 
able ordinary differentia l equation solver SUNDIALS CVODE 
(jHindmarsh et alJ I2OO5I 1 to compute the temperature and 
abundances for a set of non-linearly coupled chemical and 
thermal rate equations. The evolved species consist of H, H^, 
H~, H+, H2, He, He+, H e ++, D, D+, HP, an d free electrons 
(|Glover Jappsenl [ioQTl : IClark et al.1 l2011al l. The thermal 
energy of the gas is evolved alongside the other rate equa- 
tions. The most important coolant at low densities is molec- 
ular hydrogen, which is collisionally excited and decays ra- 
diatively via rovibrational transitions. A second, less impor- 
tant coolant is hydrogen deuteride (HD), which may become 
relevant at low temperatures due to its permanent dipole 
moment dFlower et al.ll2000l: iN akamura & Umemura 2002; 



Grei 



2007 



Bromml I2OO6I: I Johnson Bromm 2006 : Ripamontl 
McGreer &: Brvanll2008l V In fact, two of the minihaloes 



investigated here cool to ~ 100 K through HD line emis- 
sion before becoming Jeans-unstable (see also Greif et"al] 
l201lh . As will be discussed in Section 3, their evolution dif- 
fers somewhat from the haloes that have only cooled via H2. 
We note that we neglect any reactions involving deuterium 
species above a density of 10^°cm~^, since they no longer 
affect the chemical evolution of the gas. 

Above densities of 10^ cm~^, three-body reactions con- 
vert the remaining atomic hydrogen into a fully molecular 
gas. We derive the three-body formation rate by applying 
the principle of detailed bala nce to the coll isional dissocia- 
tion rate of H2 quoted in Palla et al.l (| 19831 ) . The resulting 
rate is in termediate in terms o f the large uncertainty dis- 
cussed in 'Glove r Abel (|2008h . which is reflected by the 
substantial variation of the thermal and morphological evo- 
lution of primordial gas clouds found in Turk e t al.l (|2011|). 
For riH ^ lO^cm"^, column densities are high enough 
that the gas becomes optically thick to molecular hydro- 
gen lines, and we determine the escape fration using the 
Sobolev approximation. At densities > lO^^cm"^, collision- 
induced emission (CIE) becomes i mportant and provides 
the last radiative coolin g channel (|Omukai &: Nishil Il998l : 
iRipamonti Abellliopj ). However, i t is quickly suppressed 
by the continuum opacit y of the gas (Ripam onti et al. I I2OO2I : 
IRipamonti Abelll2004l V Finally, collisional dissociation of 
H2 removes 4.48 eV of thermal energy per molecule from the 
gas and substantially offsets the compressional hea ting up 
to very high densities (|Omukaill20Q0l : lYoshida et al.| [2008). 

Once the density exceeds 10^^ cm~^, the timestep of the 
non-equilibrium solver becomes prohibitively small, which is 
due to the dependence of the three-body reactions on the 
cube of the density. We therefore switch to an equilibrium 
solver on a cell- by-cell basis for nn > ^thresh = lO^^cm"^. 
We determine the H2 fraction according to the ratio of 
the three-body formation and collisional destruction rates, 
which in our case is given by: 



riHo 



1.05 X 10" 



- exp 



5.2 X 10^ K 



(2) 



(nH-2nH2)2 ^^^^ \^ T 

where is the number density of hydrogen molecules and 




0.1 1.0 10.0 0.1 1.0 10.0 
r [AU] r [AU] 

Figure 2. From top left to bottom right: surface density, sound 
speed, rotation rate, and Toomre Q parameter versus distance to 
the centre of the primary protostar. The various line styles denote 
the minihaloes according to the legend. The above quantities are 
calculated just before the disc fragments. The growth of density 
perturbations requires a Toomre Q parameter below unity. This 
is the case below about 0.1 AU and around 1 AU, which leads to 
the formation of spiral arms in the disc. 



T the temperature of the ga s. The H^ fraction is computed 
in a similar fashion (see also lYoshida et al]|2008l l: 



/ 2Timek^T 



hi 



I exp ( - 



Xh \ 



(3) 



where is the number density of ionised hydrogen atoms, 
rrie the electron mass, /cb Boltzmann's constant, /ip Planck's 
constant, and xh the ionization energy of hydrogen. For sim- 
plicity, we assume that the He^ abundance is equal to the 
H^ abundance, multiplied by the ratio of helium to hydro- 
gen nuclei, and set the He^^ abundance to zero. We do not 
expect our results to be sensitive to this simplification. We 
converge on the self-consistent temperature and H2 fraction 
via iteration, which yields the chemical heating and cooling 
rate of the gas by H2 formation and dissociation. 

We note that we have neglected a number of additional 
cooling mechanisms that may become important at the high 
densities encountered in our simulations. Among these are 
H~ formation cooling, Lya cooling, and collisional ionization 
cooling. A more sophisticated treatment would also require 
full radiative transfer of the various line and continuum pro- 
cesses, instead of assuming a local escape fraction. An up- 
per limit on the heating of the gas by accretio n luminosity 
from the central protostar was investigated in IClark et al.l 
(|2011bh . [&eif et al.1 (|201l'), and 'Smit h et~all (|2011a"), where 
the gas was assumed to be optically thin to blackbody emis- 
sion from the protostar. In this case the fragmentation of 
the gas in the inner few AU was suppressed. However, since 
we follow the collapse of the gas to very high densities, we 
refrain from using an optically thin treatment and instead 
use escape fractions. 

Finally, we note that the chemistry solver employed 
here assumes that the electron fraction transitions into equi- 
librium at the same density as H2, instead of at nn ^ 
lO^^cm"^, where the transition is expected to occur. We 
therefore underestimate the electron fraction once the gas 
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0.1 1.0 10.0 0.1 1.0 10.0 0.1 1.0 10.0 0.1 1.0 10.0 

r [AU] r [AU] r [AU] 



Figure 3. From top left to bottom right: density of hydrogen 
nuclei, temperature, H2 fraction, and ratio of cooling time to the 
free-fall time versus distance to the centre of the primary proto- 
star. The various line styles denote the minihaloes according to 
the legend. The above quantities are calculated just before the 
disc fragments. The temperature of the gas rises only moderately 
as long as sufficient H2 is present and dissociation cooling is ac- 
tive. Below about 1 AU, cooling becomes very inefficient and the 
ratio of the cooling time to the free-fall time increases dramati- 
cally. Fragmentation occurs on a scale of about 1 AU, where both 
the Q parameter and the time-scale ratio drop below unity. 



Figure 4. Heating and cooling processes that operate in the disc 
just before it fragments. The various line styles denote compres- 
sional heating (red solid line), expansion cooling (blue solid line), 
H2 formation heating (red dotted line), H2 collisional dissociation 
cooling (blue dotted line), collision-induced emission (blue dashed 
line), and H2 line cooling (blue dot-dashed line). Each panel shows 
a different minihalo. The most important heating process is com- 
pressional heating, which is counteracted by expansion cooling 
and H2 collisional dissociation cooling. Collision-induced emis- 
sion peaks on a scale of about 1 AU, while H2 line cooling ceases 
to be important below about 10 AU. 



exceeds the threshold density for the first time, but its abun- 
dance remains so low that it does not influence the initial 
collapse of the gas. This is no longer true once a significantly 
ionised parcel of gas expands and crosses the threshold den- 
sity in the reverse direction, which leads to an unphysical 
jump in the electron fraction. We have verified that only 
a small fraction of the mass and volume of the gas at the 
outskirts of the cloud are subject to this problem, so that 
it hardly affects our results, but it should be addressed in 
future work. 



2.6 Protostar finder 

To be able to identify protostars, we need a method that 
decides whether gas clumps are protostars, and determines 
their properties such as mass and radius over time. We 
achieve this with a post-processing algorithm that uses the 
snapshots saved by AREPO and are written every ^ 0.02 yr. 
Starting from the first snapshot, we locate the densest cell 
in the simulation with at least nu = 10^^ cm""^. This cell 
is defined as the centre of a new protostar, provided that 
the distance to the centre of any other protostar exceeds 
the radius of the protostar as well as a predefined merger 
radius rmerge, which we set to 0.1 AU. The above density 
threshold gives a good indication of when the collapse of 
the gas stalls, a shock forms, and a new protostar is created, 
and the above merger radius yields the approximate point 
in time when two clumps have moved so close together that 
a merger is inevitable. 

Once a candidate cell has been selected, we determine 
the optical depth of A^ang — 10^ nearly uniformly spaced 
angular bins and A^rad = 200 logarithmically spaced radial 
bins between 0.01 and 10 AU centred on the candidate cell: 



ATj,fc = pj,ki^j,kArk, (4) 

where j denotes the angular bin, k the radial bin, pj^k 
the mass-weighted density of the bin, Kj^k the Rosseland 
mean opacity, and Ark the radial extent of the bin. For pri- 
mordial gas in chemical and thermal equilibrium, which is 
a good approximation since the density typically exceeds 
nu = 10^^ cm~^ in the protostars, the opacity is only a 
function of density and temper ature, and has been tab- 
ulated bv iMaver &; Duschl (|2005 V We linearly interpolate 
from their tables to determine the mean opacity of each bin 
from the mass- weighted density and temperature of all cells 
contributing to the bin. 

The integrated optical depth is obtained from equa- 
tion (4) by summing up the differential optical depths from 
large to small radii: 

i^k 

^j,k = ^ Ar^- (5) 

We then find the radial index /ccrit at which the spherically 
averaged escape fraction, given by: 

o _ 1 1 -exp(-r^-fc) . . 

Pesck - ^ 2^ — (6) 

j 

drops to /3crit = 1 — exp ( — 1) — 0.63, which corresponds to 
an optical depth of unity. The photospheric radius of the 
protostar is then simply identified with rfc(A:crit), and 
the mass of the protostar M* is given by the mass enclosed 
within i?*. Note that this definition implies that individual 
protostars may overlap and contain the same cells, which 
typically results in similar protostellar masses just before a 
merger occurs. 

Once a protostar has been created, it is assigned a 



© 0000 RAS, MNRAS 000, 000-000 



Primordial protostellar systems 7 




log % [cm" ] 



Side Length: 10 AU 



12 14 



16 



1^ 



20 



Figure 5. Density projections in a cube of side length 10 AU that show the evolution of the protostellar system. Each row corresponds 
to a different minihalo. The time after the formation of the primary protostar increases from left to right. The density of hydrogen 
nuclei is weighted by the density squared along the line of sight, which lies perpendicular to the plane of the disc. The disc around the 
primary protostar fragments into a number of secondary protostars, most of which migrate towards the centre of the cloud. However, 
some also obtain angular momentum from other protostars during close encounters and migrate to higher orbits. An example is the 
leftmost protostar at the last output time in MH3. A more detailed analysis of this figure is presented in Section 3.3. 



unique ID, and in addition we save the ID's of all cells that 
lie within the merger radius of the protostar. The position 
of the protostar in the next timestep is then determined 
from the new centre of mass of the saved cells. If this posi- 
tion lies within the merger radius of another protostar, it is 
merged with that protostar. The sequence of ID's therefore 
also tracks the formation history of the protostars. After 
the merging is complete, we continue with the first step and 
again find new protostars. This procedure is repeated until 
no candidate cells remain and all protostars in the snapshot 
have been identified. Finally, we loop over all snapshots until 
the entire history of the protostellar system has been recon- 
structed. 



3 RESULTS 

3.1 Collapse of central gas cloud 

Many aspects of the collapse of primordial gas in minihaloes 
have been discussed in previous work f e.g.jAbel et al.1 l2002l : 
iBromm et al .11 20021 : IVoshida et al.l l2003. 2006, 2008). We ex- 
pand on these studies by comparing the baryonic proper- 



ties of the minihaloes investigated here over an unprece- 
dented range in scales. In Fig. 1, we show a number of phys- 
ical quantities as a function of the distance to the densest 
point in each minihalo. The haloes are denoted by differ- 
ent line styles, which have been constructed using data from 
two snapshots: the inner 10^ AU are extracted from the fi- 
nal snapshots of the initial resimulations, when the density 
of hydrogen nuclei first exceeds lO^^cm"^, and the outer 



10^ 



10' AU are taken from the final snapshot of the orig- 



inal cosmological simulations, when the density of hydro- 
gen nuclei first exceeds 10^ cm"*^. The maximum density of 
10^^ cm""^ gives a good estimate of when the collapse of the 
gas stalls and a shock forms at the protostellar surface. Due 
to the self-similar nature of the collapse, moving from large 
to small radii may also be viewed as a progression in time. 
Fig. 1 therefore also shows the evolution of the central gas 
cloud from the assembly of the halo until the formation of 
the first protostar. 

From top left to bottom right, the panels in Fig. 1 show 
the density of hydrogen nuclei, enclosed gas mass, tempera- 
ture, H2 fraction, e~ fraction, radial velocity, sound speed, 
radial Mach number, rotation velocity, ratio of rotation to 
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Figure 6. Temperature projections in a cube of side length 10 AU showing the evolution of the protostellar system. Each row corresponds 
to a different minihalo. The time after the formation of the primary protostar increases from left to right. The temperature is weighted 
by the density squared along the line of sight, which lies perpendicular to the plane of the disc. The gas that accretes onto protostars 
is shock-heated to about 10^ K. From the extent of these regions it is evident that the primary protostar is generally much larger and 
more massive than the secondary protostars. The shock-heating of the gas ahead of spiral arms and the expansion cooling behind spiral 
arms are also apparent. 



Keplerian velocity, angle between the vector of the rotation 
velocity of a given radial bin and the innermost radial bin, 
and turbulent Mach number. With a few exceptions, these 
quantities are determined by using the mass- weighted aver- 
age of the cells contributing to each bin. In panels involving 
the radial or rotation velocity, we first subtract the bulk ve- 
locity of the halo from each cell. We determine the rotation 
velocity by adding up the angular momentum of all con- 
tributing cells, and divide it by the mass and radius of the 
bin. The angle (j) is given by the angle between the vector of 
the rotation velocity of a given radial bin and the innermost 
bin. We define the turbulent Mach number as: 

^turbCs = X] ^ (^^ - '^rad - Vrotf , (7) 

i 

where Cs denotes the sound speed of the radial bin, M the 
total mass, t;rad the vector of the radial velocity, Vrot the vec- 
tor of the rotation velocity, i the index of a cell contributing 
to the bin, rrii its mass, and Vi its velocity. 

To first order, the collapse of the gas follows a Larson- 
Penston solution for an isothermal, self- gravitating gas 
cloud. However, since the temperature of the gas varies and 



the cloud possesses some degree of turbulent and rotational 
support, the radial density profiles deviate somewhat from 
the idealised profile. The scatter in densities at a given 
radius can exceed an order of magnitude, showing that a sig- 
nificant amount of substructure is present. The virialisation 
heating of the gas as it falls into the DM halo can be seen 
at about 10^ AU 50 pc) in the temperature profile. As 
the density and temperature rise, the molecular hydrogen 
fraction builds up to about 10""^ and the gas cools to about 
200 K on a scale of ^ 10^ AU 0.5 pc). During this initial 
free-fall phase, the turbulent Mach number typically exceeds 
unity. When the H2 level populations transition into equilib- 
rium and the cooling rate no longer scales with the density 
squared, the gas begins to heat up rapidly to about 1000 K, 
and the turbulent Mach number decreases again. Note that 
in simulations MHl and MH2, HD cooling becomes relevant 
and the gas cools to about 100 K. Although the tempera- 
ture offset has nearly disappeared by the time the gas has 
collapsed to about 10 AU, the radial velocity in these haloes 
is about a factor of two lower than in the other haloes, and 
they have also built up less turbulence. 

On a scale of about 10^ AU, three-body reactions 
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Figure 7. Left panel: Protostellar mass as a function of time after the formation of the primary protostar. Each line denotes an individual 
protostar, and the thick dotted lines denote the total mass of the protostellar system. Mergers are characterised by the equalisation of two 
lines before one disappears. This aspect of the merging process was not captured by studies that employed sink particles, since protostars 
tend to spatially overlap for some time before they merge. The figure shows that merging is highly efficient and typically occurs between 
secondary protostars and the primary protostar. This facilitates the growth of the primary protostar to typically five times the mass of 
the second most massive protostar, and limits the number of secondary protostars that survive. Right panel: Protostellar radius versus 
time. The radii of the protostars are in the expected range of 10 — 200 R©, with the primary protostar typically at least twice as large 
as the secondary protostars. 
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rapidly convert the atomic hydrogen into a fully molecular 
gas. The resulting heating is relatively mild, since the grow- 
ing H2 abundance also ensures additional cooling. Once the 
gas becomes optically thick to H2 cooling on a scale of about 
100 AU, the temperature profile steepens again, and the gas 
heats up to a few thousand Kelvin. The inner ~ 1 AU then 
begin to cool via CIE and H2 collisional dissociation, which 
mitigates the rise in temperature as the gas further con- 
tracts. After most of the H2 has been destroyed, the gas can 
no longer cool and heats up nearly adiabatically to almost 
10^ K, resulting in a rapid increase of the electron fraction. 

The rotation velocity remains comparable to the infall 
velocity throughout the entire evolution of the gas cloud, 
and its ratio to the Keplerian velocity varies between about 
0.2 and 0.7, indicating a high degree of rotational support. 
The peak of the profile is at 10^ AU after a period of 
strong heating, when the ratio of cooling time to dynami- 
cal time is large and the gas has time to settle into a disc. 
The rotation angle varies significantly during the collapse of 
the gas, which indicates vigorous transport of angular mo- 
mentum. The main agent for angular momentum transport 
during the initial collapse is turbulence, which is comparable 
in magnitude to the Mach number of the radial infall, and 
typically of order unity. 



3.2 Disc formation and fragmentation 

After the first protostar has formed, the gas becomes fully 
rotationally supported in a Keplerian disc. The stability of 
the disc to perturbations in th e density can be quantified 
with the Toomre Q parameter (|Toomrelil963 ): 



where n is the epicyclic frequency of the disc, S the surface 
density, and Cs the sound speed of the gas. We replace the 
epicyclic frequency n with the orbital frequency Q, which 
is appropriate for Keplerian discs. The Q parameter deter- 
mines whether density perturbations in an infinitely thin, 
isothermal disc can grow. For thick discs, as is the case here, 
a similar criterion may be found that devi ates only by a fac- 
tor of order unity (e.g. IWang et al.ll201Ql l. For Q > 1, the 
pressure of the gas and the shear by the differential rotation 
of the disc are sufficient to prevent local collapse, while for 
Q < 1 the disc fulfils the Toomre criterion and becomes un- 
stable, leading to the formation of spiral arms that transport 
mass inwards and angular momentum outwards. 

In Fig. 2, we show the surface density, sound speed, or- 
bital frequency, and Q parameter in mass- weighted spherical 
shells around the densest cell in each minihalo just before 
the first fragment forms. We note that applying the same 
analysis exclusively to cells in the disc gives very similar re- 
sults, since the mass within radial shells is dominated by the 
disc component. In the central few AU, the surface density is 
much higher as compared to studies that have e mployed sink 
particl es, and does not show a central dip (e.g. Clar k et al.l 
l2011br ). This is not surprising, since these studies typically 
lose resolution on the scale of the accretion radius. The Q 
parameter therefore does not diverge at small radii, but de- 
creases to well below unity, which is reflected by the devel- 
opment of pronounced spiral arm patterns on sub-AU scales 
at very early times. The orbital frequency on these scales 



- MHl 

MH2 




2 4 6 8 10 
t[yr] 

Figure 8. Total accretion rate of the protostellar system versus 
time after the formation of the primary protostar. The various 
line styles denote the minihaloes according to the legend. The 
accretion rates vary strongly between about 10~^ and 1 yr"-*^, 
with a mean of about O.lM0yr-i. In MHl and MH2, HD cooling 
becomes important and leads to a systematically lower accretion 
rate. 

drops to well below 10yr~^. The rapid decline of the sound 
speed leads to a second minimum at ^ 1 AU, followed by an 
increase at larger radii due to the flat surface density profile 
and nearly constant sound speed. Interestingly, the gas does 
not fragment at the location of the first minimum, but at 
the second minimum at about 1 AU. This shows that Q ^ 1 
is not a sufficient criterion for gravitational fragmentation. 

A commonly employed criterion in addition to the 
Toomre criterion to e valuate whethe r the gas fragments is 
the Gammie criterion (|Gammiell200lh : 
3 

^'cool < T^torb, (9) 

where tcooi is the cooling time and torb is the orbital time. 
Again, this equation strictly applies only to infinitely thin 
discs, but a similar criterion may be found for thick discs. In 
Fig. 3, we show the density of hydrogen nuclei, temperature, 
H2 fraction, and ratio of cooling time to free-fall time using 
the same output times and binning as in Fig. 2. Note that 
the orbital time may be approximated by the free-fall time, 
since they differ only by a factor of order unity in a Keplerian 
disc. Moving from large to small radii, the ratio of cooling 
time to free-fall time rises sharply above unity once most 
of the H2 has been dissociated and the gas can no longer 
cool, which is accompanied by a rapid rise in temperature 
as the gas heats up to > 10^ K. This occurs on a scale of 
about 1 AU, which corresponds to the second minimum in 
the Q parameter. This is the characteristic radius at which 
the disc fragments. 

In Fig. 4, we investigate the various heating and cooling 
processes that result in the nearly isothermal radial profile 
of the disc down to very small radii. The rates are computed 
per unit mass instead of unit volume due to the strong in- 
crease in density towards the centre. We use the same out- 
put times and binning as in Figs. 2 and 3. The individual 
lines correspond to compressional heating, expansion cool- 
ing, three-body H2 formation heating, H2 collisional disso- 
ciation cooling, CIE cooling, and H2 line cooling. The most 
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Figure 9. From top left to bottom right within each panel: density of hydrogen nuclei, temperature, sound speed, radial velocity, rotation 
velocity, ratio of radial to rotation velocity, ratio of rotation to Keplerian velocity, and turbulent Mach number versus distance to the 
centre of the primary protostar. Each panel shows an individual minihalo, and the various lines styles denote the above quantities at the 
output times shown in the legend. The grey lines in the top left panels indicate the radii of the primary protostars in each halo. This 
figure is discussed in detail in Section 3.3. 
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Figure 10. Average number of secondary protostars formed 
(solid line), merged (dotted line), and merged with the primary 
protostar (dashed line) in each halo. At the last common out- 
put time, about half of the secondary protostars have merged 
with the primary protostar, and only about one third survive. 
Nevertheless, the multiplicity of the protostellar system increases 
monotonically until the end of the simulation. 

important heating process is compressional heating, which 
typically dominates over H2 formation heating by an order 
of magnitude. The cooling of the gas is dominated by expan- 
sion cooling and collisional dissociation of H2, which nearly 
balance the compressional heating rate. The CIE cooling 
curve peaks on a scale of about 1 AU, but otherwise has a 
relatively small contribution to the total cooling rate. H2 
line cooling is important on scales ^10 AU, but is quickly 
diminished by the opacity of the gas on smaller scales. In 
MH3 and MH4, the dissociation cooling rate drops sharply 
at 0.1 AU, since the molecular hydrogen reservoir is ex- 
hausted on larger scales as compared to MHl and MH2 (see 
Fig. 3). 

The above analysis shows that even though perturba- 
tions in the disc grow on scales as small as < 0.1 AU, the 
disc does not fragment on these scales since the gas is opti- 
cally thick and the cooling time is larger than the dynami- 
cal time. Instead, fragmentation occurs on a scale of about 
1 AU, where H2 collisional dissociation and collision-induced 
emission are sufficient to reduce the cooling time below the 
dynamical time. This is nearly an order of mag nitude below 
the sca le on which the first fragments formed in I Clark et al.l 
(l2011bl l. where the gas was assumed to be optically thin 
to radiation from the primary protos tar and fragmenta tion 
was suppressed below about 10 AU. In Greif et al.l (|201 lh. an 
indication for a reduced fragmentation scale was presented 
in section 3.5, where in addition to the reduced fragment 
masses shown in fig. 16, the distance of the first fragment to 
the primary decreased as the accretion radius was reduced. 



3.3 Evolution of protostellar system 

The fragmentation of the disc into one or more secondary 
protostars further complicates the evolution of the central 
gas cloud. This is evident from Fig. 5, where we show the 
number density of hydrogen nuclei in cubes with a side 
length of 10 AU at six different output times. The line of 



sight is aligned perpendicular to the plane of the disc. The 
initial burst of star formation after the disc becomes grav- 
itationally unstable is followed by a complicated succession 
of mergers between existing protostars and the formation of 
new protostars. By the end of the simulation, a small system 
of protostars ranging from a binary in MHl to a system of 
five protostars in MH3 and MH4 have formed. 

The temperature of the gas is shown in Fig. 6. The 
hottest parcels of gas reside within protostars and have tem- 
peratures in excess of 10^ K. The more massive protostars, 
and in particular the primary protostar, tend to reside at the 
centre of the cloud and are hotter than their low-mass coun- 
terparts at the outskirts of the cloud. The compressional 
heating of the gas ahead of spiral arms and the expansion 
cooling behind spiral arms are also apparent. Some striking 
differences between haloes that have activated HD cooling 
and those that have not are evident. For example, the gas 
in simulations MHl and MH2 is significantly less turbulent, 
which is reflected by the s moother appearance of the density 
and temperature (see also Clark et al.l l2011al V The discs in 
these haloes remain stable for a longer period of time, and 
the gas tends to be more condensed, which is evident from 
the smaller sizes of the discs and the protostars. 

The evolution of the masses and radii of the protostars 
are shown in Fig. 7. Each line denotes an individual proto- 
star, and the total mass of all protostars is denoted by the 
thick dotted line. Including the primary protostar, a total 
of 5, 10, 9, and 11 protostars form, out of which 2, 4, 5, 
and 5 survive until the end of the simulation. The masses 
of the protostars range from a few times 10~^ to almost 
IM0. The total mass is dominated by the primary proto- 
star at the centre of the cloud, which is typically five times 
as massive as the second most massive protostar. The radii 
of the protostars are shown in the right panel of Fig. 7, and 
are in th e expected range of a few tens to two hundred so- 
lar radii teosokawa Omukaillioool : ISmith et al.ll2011bl ). In 
Fig. 8, we show the total accretion rate of all protostars, av- 
eraged over a one- year time period. The accretion rate varies 
strongly between about 10~^ and IMoyr"^, with a mean 
of about 0.1 M© yr~\ 

In Fig. 9, we analyse the central gas cloud as a whole 
as the protostellar system evolves. From top left to bot- 
tom right, the panels show the density of hydrogen nuclei, 
temperature, sound speed, radial velocity, rotation velocity, 
ratio of radial to rotation velocity, ratio of rotation to Kep- 
lerian velocity, and turbulent Mach number as a function of 
distance to the centre of the primary protostar. The various 
line styles denote the output times according to the legends, 
and the grey lines also show the radii of the primary pro- 
tostars. The profiles are centred on the first protostar and 
have been determined according to the description in Sec- 
tion 3.1, but using only data from the final resimulations. 
The spatial range is also much smaller, since only the cen- 
tral regions have had enough time to evolve from their initial 
state. The formation of additional protostars over time is ev- 
ident from the bumps in the density profiles, which show an 
increasing deviation from the initial profiles as the gas clouds 
become highly inhomogeneous. The temperature beyond the 
primary protostar does not increase by much, showing again 
that cooling within the disc is very efficient. The location of 
the upturn of the radial velocity, as well as the downturn of 
the ratio of the radial velocity to the rotation velocity in- 
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Figure 11. The evolution of the two most long-hved secondary protostars in each minihalo, which in addition merge with the primary 
protostar before the end of the simulation. The x-axis denotes the ratio of the time after their formation to the time that they survive. 
From top left to bottom right, the sub-panels show the distance to the centre of the primary protostar, the ratio of the angular momentum 
per unit mass to the initial value when the protostars were formed, and the ratio of the time-scales for gravitational and pressure gradient 
torques to operate to the times that the protostars survive. The thick grey dashed lines denote a ratio of unity, and the red diamonds 
and blue crosses denote torques that increase and decrease the angular momentum of the protostars, respectively. The predominance 
of crosses in the bottom left panels shows that gravitational torques tend to decrease the angular momentum of the protostars. On 
the other hand, torques exerted by pressure gradients typically increase the angular momentum of the protostars. The time-scales for 
gravitational torques to operate are generally close to the time the protostars survive, and significantly smaller than the time-scales for 
pressure gradient torques to operate, showing that gravitational torques are indeed responsible for the migration of the protostars to the 
centre of the cloud. 



creases over time and indicates where the infall stalls and the 
gas becomes rotationally supported. The Keplerian velocity 
within this region increases to nearly unity and shows the 
formation and growth of a disc that encompasses the sec- 
ondary protostars with their own small discs. The turbulent 
Mach number of the gas within this disc increases by about 
a factor of two, and is generated by the self- gravity of the 
gas and the non-axisymmetric nature of the system. Below 
about 1 AU, which is typically located within the protostars, 
the turbulent motions become subsonic. 



3.4 Merger and survival rate 

One of the key parameters that controls the stellar mass 
function is how efficiently the protostars merge with each 
other. In Fig. 10, we compare the average number of sec- 
ondary protostars formed per halo with the average number 
of mergers. By the end of the simulation, nearly two thirds 
of the secondary protostars have merged away, and only one 
third survive. The dashed line shows that about half of the 
secondary protostars have merged with the primary proto- 
star. Evidently, they shed their angular momentum very ef- 
ficiently, for which either gravitational or pressure gradient 
torques are responsible. In the following, we investigate in 
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Figure 12. The evolution of the two most long-hved secondary protostars in each minihalo, which in addition survive until the end of 
the simulation (note that in MHl only one secondary protostar survives). The x-axis denotes the ratio of the time after their formation 
to the time that they survive. From top left to bottom right, the panels show the distance to the centre of the primary protostar, the 
ratio of the angular momentum per unit mass to the initial value when the protostars were formed, and the ratio of the time-scales for 
gravitational and pressure gradient torques to operate to the times that the protostars survive. The thick grey dashed lines denote a 
ratio of unity, and the red diamonds and blue crosses denote torques that increase and decrease the angular momentum of the protostars, 
respectively. In comparison with Fig. 11, diamonds tend to dominate, which indicates that gravitational torques increase the angular 
momentum of some of the protostars, and often results in their migration to higher orbits. 



turn the evolution of secondary protostars that merge with 
the primary protostar, and of secondary protostars that sur- 
vive until the end of the simulation. 



In our first attempt to determine the various torques 
acting on the protostars, we assumed that they may be ap- 
proximated as point masses. However, this proved to be in- 
accurate, since the complex density and temperature profiles 
near the protostellar surface strongly affect the torques. This 
prompted us to use AREPO to generate an additional output 
for each snapshot that contains the gravitational accelera- 
tion and pressure gradient of all cells. We then determine the 
total gravitational and pressure gradient torques per unit 
mass Tgrav and Tpres on the protostars via summation over 
all cells that lie within the protostars: 



Tgrav = ^r'^X(m^a^), (10) 

i 

where % is the index of the cell, ri its distance to the centre 
of the primary protostar, rrii its mass, pi its density, ai its 
gravitational acceleration, and {VP)i its pressure gradient. 
We also determine the angular momentum per unit mass /, 
which is given by: 

^^r,x{m,v^), (12) 

i 

where Vi denotes the velocity of the cell with respect to 
the velocity of the primary protostar. From the angular mo- 
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Figure 13. Average survival time (solid line), free-fall time (dot- 
ted line), and dynamical friction time (dashed line) of the proto- 
stars that have merged with the primary protostar before the end 
of the simulation, shown as a function of distance to the centre of 
the primary protostar. The redistribution of angular momentum 
by gravitational torques is so efficient that it leads to the migra- 
tion of about half of the secondary protostars to the centre of the 
cloud in only a free-fall time, which is the smallest time-scale on 
which a gravitationally driven process may operate. The dashed 
line shows that the torques are at least in part due to dynamical 
friction - especially at small radii where the densities are high. 



menta and torques, we determine the time-scales for gravi- 
tational and pressure gradient torques to operate: 

|2 



^erav — 



\l\ 



and 



\i\ 



(13) 



Note that positive values indicate that the angular momen- 
tum per unit mass is increasing, while negative values indi- 
cate that it is decreasing. 

In Fig. 11, we analyse the evolution of the two most 
long-lived secondary protostars in each halo, which in ad- 
dition merge with the primary protostar. The x-axis shows 
the ratio of the time after their formation to the time they 
survive. The top left panels show the distance of the proto- 
stars to the centre of the primary protostar. The top right 
panels show the ratio of the angular momentum per unit 
mass to the initial value when the protostar is formed. Fi- 
nally, the bottom panels show the ratio of the time-scales 
for gravitational and pressure gradient torques to operate 
to the times the protostars survive. Red diamonds denote a 
positive value, and blue crosses denote a negative value. For 
better visibility, we have smoothed the lines in the bottom 
two panels over a small period of time. 

A number of important conclusions may be drawn from 
this figure. First, the bottom panels show that gravitational 
torques are more important than pressure gradient torques, 
since the respective time-scales are significantly smaller. Sec- 
ond, the predominance of negative values in the bottom left 
panels shows that gravitational torques typically lead to a 
decrease of the angular momentum of the protostars, while 
the opposite is the case for the torques exerted by pressure 
gradients. Finally, the time-scales for gravitational torques 
to operate are generally close to the survival time, showing 
that gravitational torques are indeed responsible for the mi- 
gration of the protostars to the centre of the cloud. A closer 



look at the evolution of individual protostars also shows that 
the dips and peaks in the time-scale profiles correlate well 
with those in the distance and angular momentum. For ex- 
ample, the solid line in MH2 shows a clear peak in the dis- 
tance and angular momentum at about t/tsurv = 0.5, which 
is accompanied by a momentary increase of tgrav to positive 
values. In a second example, the distance and angular mo- 
mentum of the protostar denoted by the dotted line in MH3 
continuously decline, which is reflected by an extended pe- 
riod of negative gravitational torques. Although not entirely 
apparent from the bottom panels, due to the smoothing ap- 
plied to the profiles, the torques acting on the protostars are 
also highly variable in time. 

In Fig. 12, we show the evolution of the two most long- 
lived secondary protostars in each halo, which in addition 
survive until the end of the simulation (note that in MHl 
only one secondary protostar survives). In this case, and in 
contrast to Fig. 11, gravitational torques tend to be domi- 
nated by positive values, and may lead to an increase of the 
distance and angular momentum of individual protostars. 
Two examples are the protostars denoted by the dotted line 
in MH2 and the solid line in MH4, which increase their sep- 
aration to the centre of the primary protostar by a factor 
of a few. This is caused by torques from nearby protostars, 
indicated for example by the correlation of the dotted line 
in MH2 in Fig. 12 with the solid line in Fig. 11. The result- 
ing migration of the protostars to higher orbits is si milar 
to the dynamical ejections found in iGreif et al ] (|201lh . but 
less severe since a fraction of the angular momentum of the 
protostars is transferred to the surrounding gas instead of 
to other protostars. 

In Fig. 13, we compare the average free-fall time and dy- 
namical friction time of the secondary protostars that merge 
with the primary protostar to the time they survive. The 
free-fall time tff is given by: 



tff 



Stt 



32Gpc 



(14) 



where G is the gravitational constant and pcioud is the mean 
density of the gas, which we estimate by using the mass 
enclosed within the distance of the protostar to the centre 
of the primary protostar. We also determine the dynamical 
friction time tfric, which is approximately given by: 



tfri 



(15) 



where Ffric is the dynamical friction force. We estimate the 
dynamical friction force with the approximate expression for 
collisionless systems: 



47tG^M^PcIo 



(16) 



where is the velocity of the protostar with respect to 
the velocity of the primary protost ar. We note th at more 
sophisticated expressions exist fe.g. IOstrikerlll999h . but in 
light of the complicated morphology of the gas cloud their 
applicability is questionable here. The above time-scales are 
calculated using data from all haloes and output times, and 
we bin them according to the distance of the protostars to 
the centre of the primary protostar. 

Fig. 13 shows that the reorganization of the angular mo- 
mentum of the protostars and the resulting migration to the 
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centre of the cloud typically occur in a free-fall time, which 
is the smallest time-scale on which a gravitationally driven 
process may operate. As indicated by the dashed line, dy- 
namical friction may be partially responsible for the strong 
gravitational torques, in particular at small radii where the 
density is high. Other agents include irregularities in the 
gas cloud induced by turbulence or spiral arm patterns, and 
nearby protostars. A more sophisticated analysis of the ori- 
gin of the torques would require an implementation within 
AREPO of the protostar finder described in Section 2.6, to- 
gether with an auxiliary routine to decompose the gravi- 
tational force into its various components. However, this is 
beyond the scope of the present study. 



4 SUMMARY AND CONCLUSIONS 

Recent simulations of primordial star formation have sug- 
gested that the gas at the centre of minihaloes fragments 
into a system of protostars, as a consequence of the high ac- 
cretion rate from the envelope onto the disc, and the efficient 
cooling of the disc by molecular hydrogen (Clark et al. 2008, 
2011a,b; Stacy et al. 2010, Greif et al. 2011). These stud- 
ies employed sink particles to avoid the Courant constraint, 
which allowed a continuation of the simulations beyond the 
formation of the first protostar, at the cost of not accurately 
resolving the hydrodynamics near the accretion radius. We 
have attempted to overcome this limitation by not insert- 
ing sink particles, and instead directly resolving individual 
protostars and their interaction with the surrounding gas. 
In analogy to previous work, we find that the Keplerian disc 
formed around the primary protostar fragments into a num- 
ber of secondary protostars, which is enabled by collisional 
dissociation cooling of H2 and collision-induced emission. 

Once the first fragments have formed, aggressive grav- 
itational torques between the protostars and the gas begin 
to redistribute their angular momentum. This leads to the 
migration of about half of the secondary protostars formed 
in the disc to the centre of the cloud in a free-fall time, where 
they merge with the primary protostar and facilitate its 
growth to about five times the mass of the second most mas- 
sive protostar. The preferential merging with the primary 
protostar is akin to the 'runaway growth' of a single object 
encountered in studies of planet formation (|Kokubo Idal 
1 19961 ) - as opposed to the ' oligarchic growth' of multiple ob- 
jects feokubo k Ida' ll998l ) - and to what is seen in some 
simulations of massi ve star fo rmation in the present-day 
Universe (e.g. Krum holz et al.l 120091 ). Next to the inward 
migration of a large fraction of the protostars, a few also 
receive angular momentum from other protostars due to 
N-body interactions and mig rate to high er orbits. Similar 
behavior was found in .Greif et al.l (|201lh . but with higher 
eject a velocities since the angular momentum of the proto- 
stars was transferred solely to the merger product. Here, the 
complicated torques acting on the extended protostars re- 
sults in the transfer of a fraction of the angular momentum 
to the surrounding gas. 

The relatively high merger rate leads to the survival 
of only every third protostar formed in the disc. Neverthe- 
less, the number of protostars present at any given time in- 
creases monotonically until the end of the simulation, sug- 
gesting that the protostellar system will continue to grow 



beyond the limited period of time simulated here. Despite 
the considerable computational effort involved, we have fol- 
lowed only a small fraction of the total accretion time, so 
that we cannot predict the final mass function of the first 
stars. In addition, we have not accurately solved for the ra- 
diation emitted by the protostars, which could significantly 
affect the chemical and thermal evolution of the gas (Smith 
et al. 2011a,b; Hosokawa et al. 2011; Stacy et al. 2011). An- 
other complication is given by magnetic fi elds, which are 



thought to be important in minihaloes (e.g. Machida e t al 
2008; Xu et al . 2008; Schleicher et al. 2009, 2010; Sur et al 
2OIOI : iFederrath et al.ll2011l : ISchober et al.ll2012l : ITurk et al 
2OI2I ). Next-generation simulations should therefore include 
additional physics, as well as continue the simulations to 
later times. The latter requires an increase of the raw com- 
putational speed of AREPO, which may be achieved by opti- 
mising time-consuming calculations such as the construction 
of the Voronoi tesselation, or by utilising multi-threaded en- 
vironments. In the near future, it may thus become possible 
to significantly extend the period of time that can be simu- 
lated. 
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